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Abstract 

Circulant matrices play a central role in a recently proposed formulation of three-way 
data computations. In this setting, a three-way table corresponds to a matrix where 
each "scalar" is a vector of parameters defining a circulant. This interpretation provides 
many generalizations of results from matrix or vector-space algebra. We derive the power 
and Arnoldi methods in this algebra. In the course of our derivation, we define inner 
products, norms, and other notions. These extensions are straightforward in an algebraic 
sense, but the implications are dramatically different from the standard matrix case. For 
example, a matrix of circulants has a polynomial number of eigenvalues in its dimension; 
although, these can all be represented by a carefully chosen canonical set of eigenvalues 
and vectors. These results and algorithms are closely related to standard decoupling 
techniques on block-circulant matrices using the fast Fourier transform. 

Keywords: block-circulant, circulant module, tensor, FIR matrix algebra, power 
method, Arnoldi process 



1. Introduction 

We study iterative algorithms in a circulant algebra, which is a recent proposal for a 
set of operations that generalize matrix algebra to three-way data ( Kilmer et al. . 20081 ). 



In particular, we extend this algebra with the ingredients required for iterative methods 
such as the power method and Arnoldi method, and study the behavior of these two 
algorithms. 

Given an m x n x k table of data, we view this data as an m x n matrix where each 
"scalar" is a vector of length k. We denote the space of length- k scalars as Kfe. These 
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scalars interact like circulant matrices. Circulant matrices are a commutative, closed 
class under the standard matrix operations. Indeed, is the ring of circulant matrices, 
where we identify each circulant matrix with the k parameters defining it. 

Formally, let a G K^. Elements in the circulant algebra are denoted by an underline 
to distinguish them from regular scalars. When an element is written with an explicit 
parameter set, it is denoted by braces, for example 



a = { c 



4- 



In what follows, we will use the notation -R- to provide an equivalent matrix-based no- 
tation for an operation involving Wc define the operation circ(-) as the "circulant 
matrix representation" of a scalar: 



a -H- circ(a) 



CH2 cti 



«2 



a 2 ol\ 



(1) 



Let a be as above, and also let /JeKj. The basic addition and multiplication operations 
between scalars arc then 



a + (3 <H> circ(a) + circ(/3) and a o f3 <-» circ(a) circ(/3). 



(2) 



We use here a special symbol, the o operation, to denote the product between these 
scalars, highlighting the difference from the standard matrix product. Note that the 
element 

1 = {l o ... o} 

is the multiplicative identity. 

Operations between vectors and matrices have similar, matricized, expressions. We 
use KjJ to denote the space of length-n vectors where each component is a fc-vector in Kfc, 
and Kr x " to denote the space of to x n matrices of these k- vectors. Thus, we identify 



each m x n x k table with an element of K^ xn . 



Let A e K™ x ™ and x e K r ^. Their 



product is: 



A o x = 



circfAi.i) ... circfAi, 



c±Tc(A m; i) ... circ(4 m|Tl ) 



circ(a^i) 



= (2n) 



Thus, we extend the operation circ to matrices and vectors of K k scalars so that 

circ(_Aii) ... circ(_Ai i?l ) 



circ(A) 



and circ(x) 



The definition of the product can now be compactly written as 

ioxf* circ(j4) circ(x). 



circ(xi) 



(3) 



(4) 



(5) 



Of course this notation also holds for the special case of scalar- vector multiplication. Let 
a e K fe . Then 

x o a O circ(x) circ(a). 
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The above operations define the basic computational routines to tre at m x n x k 
arrays as m x n matrices of They are equivalent to those proposed bv lKilmer et al. 
( 20081 ). an d they constitute a module over vectors composed of circulants, as shown 
recently in iBramanl ( To appear ). Based on this analysis, we term the set of operations 
the circulant algebra. We note that these operations have more efficient implementations, 
which will be discussed in Sections [3] and [6] 

The circ ulant algebra an alyzed in this paper is closely related to the FIR matrix al- 
gebra due to lLambert (|l996l Chapter 3). Lambert proposes an algebra of circulants; but 
his circulants are padded with additional zeros to better approximation a finite impulse 
respo nse operator. He uses it to study blind deconvolution problems ([Lambert et al 



20011 ). As he observed, the relationship with matrices implies that many standard de- 



compositions and techniques from real or complex valued matrix algebra carry over to 
the circulant algebra. 

The circulant algebra in this manuscript is a p articular instance of a matrix-over-a - 
ring, a long studied generalization of linear algebra ( Mcdonald . 1984 ; Brewer et all 
Prior work focuses o n Roth theore ms for the equation AX — XB = C (IGustafsonl. 



1986). 



1979): 



gener alized inverses (jPrasa d, 1994); completion and controllability pr oblems dGurvits et al 



1992); matrices over the ring of integers for computer algebra systems Earner and McCurlev 



1991 ); and trans fer functions and linear dynamic systems (jSontael Il976l) . Finally, see 



Gustafsonl (|199l[ ) for some interesting relationships between vecto rs space theo r y and 



module theory. A recent proposal ext ends many of the ope rations in iKilmer et ail (|2008l) 



to more general algebraic structures (jNavasca et al. , 20101 ) 



Let us provide some further context on related work. Multi-way arrays, ten sors, 
and hypermatrices are a burgeoning area of research; see iKolda and Bader (2009) for 
a recent comprehensive survey. Some of the major themes are multi-linear ope r ations , 
fitting multi- linear models, and multi-linear generalizations of eigenvalues ([Oil . 12007 ). 
The formulation in this paper gives rise to stronger relationships with the litera ture 
on blo ck-circulant matrices, which have been studied for quite some time. See I Tee 
(|2005l ) and the references therei n for further historical and mathematical context on 
circulant matrices. In particular, I Baker dl989h gives a procedure for the SVD of a block 
circulant that involves using the fast Fourier transform to decouple the problem into 
independent sub-problems, just as we shall do throughout this manuscript. Other work 
in this vein includes solving block-circulant systems that arise in the theory of antenn a 
arrays: (jSinnott and Harrington! . 1 1973tlDe Mazancourt and Gerlidll983[IVescovcil . ll997l ). 

The remainder of this paper is structured as follows. We first derive a few necessary 
operations in Section [21 including an inner product and norm. We then continue this dis- 
cussion by studying these same operations using the Fourier transform of the underlying 
circulant matrices (Section |3|). A few theoretical properties of eigenvalues in the circulant 
algebra are analyzed in Section 21 Th at section is a necessary prelude to the su bsequent 
discussion of ho w the power method (|von Mises and Pollaczek-Geiringerl Il929l ) and the 
Arnoldi method ( Krvlovl 1 1931 ; Lanczosl Il950l: Arnoldi . 1951 ) generalize to this algebra, 
which comes in Section [5] We next explain how we implemented these operations in 
a Matlab package (Section [6]) ; and we provide a numerical example of the algorithms 
(Section [7]). Section [8] concludes the manuscript with some ideas for future work. 
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Figure 1: The power method for a matrix A 6 R 

Require: A,x.(°\t 

1: xC ><-x(°) ||x(°)||" 1 

2: for k = 1, . . . , until convergence do 
3: y( fe) 4- Ax (fc_1 ) 

4: <— lly^ll 

5: xW^yWaW" 1 

6: if ||sign(x^)x( fc ) — sign(xj fc ~^)x( fc_1 )|| < r then 

7: return x' fc ^ 

8: end if 

9: end for 



2. Operations with the power method 

In the introduction, we provided the basic set of operations in the circulant algebra 
(eqs. (H}-©). We begin this section by stating the standard power method, and then 
follow by deriving the operations it requires. 

Let A £ R raxn and let x g K™ be an arbitrary starting vector. Then the power 
method proceeds by repeated applications of A; see Figure Q] for a standard algorithmic 
description. (Line [6] checks for convergence and is one of sever al possible stopping crite- 



ria.) Under mild and well-known conditions (see lStewartl (|2001I )L this iteration converges 
to the eigenvector with the largest magnitude eigenvalue. 

Not all of the operations in Figure Q] are defined for the circulant algebra. In the first 
line, we use the norm ||x(°)|| that returns a scalar in R. We also use the scalar inverse 
a . The next operation is the sign function for a scalar. Let us define these operations, 
in order of their complexity. In the next section, we will reinterpret these operations in 
light of the relationships between the fast Fourier transform and circulant matrices. This 
will help illuminate a few additional properties of these operations and will let us state 
an ordering for elements. 

2.1. The scalar inverse 

We begin with the scalar inverse. Recall that all operations between scalars behave 
like circulant matrices. Thus, the inverse of a € is 

a^ 1 ■<-> circ(a) . 

The matrix circ(a) -1 is also circulant ( Davisl Il979). 



2.2. Scalar functions and the angle function 

Other scalar functions are also functions of a matrix (see Higham (2008)). Let / be 
a function, then 

f(a) O /(circ(a)) 

where the right hand side is the same function applied to a matrix. (Note that it is not 
the function applied to the matrix element- wise.) 
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then angle(re ) = e %b . For real or complex numbers x, we then 



The sign function for a matrix is a special case. As explained in Higha m (2008), the 
sign function applied to a complex value is the sign of the real- valued part. We wish to 
use a related concept that generalizes the real-valued sign that we term "angle." Given 
a complex value re 
have 

angle(a;) |af| = x. 

Thus, we define 

angle(a) f-> circ(abs(a)) _1 circ(a) 



2.3. Inner products, norms, and conjugates 

We now proceed to define a norm. The norm of a vector in KJJ produces a scalar in 

K fe : 

/ n \ 1 / 2 

||x|| O (circ(x)* circ(x)) 1/2 = ^ circ(gj)* circ(xi) I 

For a standard vector x £ C™, the norm ||x|| = \/x*x. This definition, in turn, follows 
from the standard inner product attached to the vector space C n . As we shall see, our 
definition has a similar interpretation. The inner product implied by our definition is 

(x, y) O circ(y)* circ(x). 

Additionally, this definition implies that that the conjugate operation in the circulant 
algebra corresponds to the transpose of the circulant matrix 

a O circ(a)*. 

With this conjugate, our inner product satisfies two of the standard properties: conjugate 
symmetry (x, y) = (y,x) and linearity («ox,y) — a o (x, y). The notion of positive 
definiteness is more intricate and we delay that discussion until after introducing a de- 
coupling technique using the fast Fourier transform in the following section. Then, in 
Section 13.31 we use positive definiteness to demonstrate a Cauchy-Schwarz inequality, 
which in turn provides a triangle inequality for the norm. 



3. Operations with the fast Fourier transform 

In Section [2j we explained the basic operations of the circulant algebra as operations 
between matrices. All of these matrices consisted of circulant blocks. In this section, we 
show how to accelerate these operations by exploiting the relationship between the fast 
Fourier transform and circulant matrices. 

Let C be a k x k circulant matrix. Then the eigenvector matrix of C is given by the 
k x k discrete Fourier transform matrix F, where 

F- - JL w (<-UCj'-i) 

V Vk 

and lo = e 2n ^ k . This matrix is complex symmetric, F T — F, and unitary, F* = F^ 1 . 
Thus, C — FDF* , D = diag(Ai, . . . , Afe). Recall that multiplying a vector by F or F* 
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circ(A) 



:(A)(I® F) 



eft (A) 



Figure 2: The sequence of transformations in our eft operation. Given a circulant A, we convert it 
into a matrix by circ(A). The color of the circles in the figure is emphasizing the circulant structure, 
and not equality between blocks. In the third figure, we diagonalize each circulant using the Fourier 
transform. The pattern of eigenvalues is represented by squares. Here, we are coloring the squares to 
show the reordering induced by the permutation at the final step of the eft operation. 



can be accomplished via the fast Fourier transform in O(fclogfc) time instead of 0(k 2 ) 
for the typical matrix-vector product algorithm. Also, computing the matrix D can be 
done in time O(fclogfc) as well. 

To express our operations, we define a new transformation, the "Circulant Fourier 
Transform" or eft. Formally, eft : a £ Kk i-» £ kxk an d its inverse icft : £ kxk i-» Kk 
as follows: 



cft(a) 



F* circ(a)F, icft 



: a Fc±t(g)F*, 



where 6tj are the eigenvalues of circ(a) as produced in the Fourier transform order. 
These transformations satisfy icf t(cf t(a)) = a and provide a convenient way of moving 
between operations in Kk to the more familiar environment of diagonal matrices in <C kxk . 

The eft and icft transformations are extended to matrices and vectors over 
differently than the circ operation we saw before. Observe that eft applied "element- 
wise" to the circ(A) matrix produces a matrix of diagonal blocks. In our extension of the 
eft routine, we perform an additional permutation to expose block-diagonal structure 
from these diagonal blocks. This permutation P m AP^ transforms an mk x nk matrix 
of k x k diagonal blocks into a block di agonal mk x nk with m x n size blocks. It is also 
known as a stride permutation matrix (jGranata et all 11992 ). The construction of P m , 
expressed in Matlab code is 

p = reshaped :m*k,k,m) ' ; 

Pm = sparse (1 :m*k,p( :), 1 ,m*k,m*k) ; 

The construction for P n is identical. In Figure [21 we illustrate the overall transformation 
process that extends eft to matrices and vectors. 

Algebraically, the eft operation for a matrix A £ K™ x ™ is 

cf t(A) = P m (I m ® F*) circ(A)(I„ ® F)P^, 



where P m and P n are the permutation matrices introduced above. We can equivalently 
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write this directly in terms of the eigenvalues of each of the circulant blocks of circQ4): 



cftQ4) 



... a: 



where Aj' s , . . . , Al' 5 are the diagonal elements of cft(.A r)S ). The inverse operation icf t, 
takes a block diagonal matrix and returns the matrix in K™ XTi : 

icft(A) o (J m <g> F)P T m AP n {I n ® F*). 

Let us close this discussion by providing a concrete example of this operation. 



Example 1. Let A = 



illustrated in Figure^ are: 



{2 3 1} {8 -2 0} 
{-2 2} {3 1 1} 



The result of the circ and eft operations, as 



circ(A) 



2 13 

3 2 1 
13 2 



8 0-2 
-2 8 
0-2 8 



-2 2 
0-2 2 
2 0-2 



3 11 
13 1 
113 



(I®F*) circ(A){I®F) = 



6 


6 






V3i 


-9-V3i 





5 


-3+V3i 


2 


-3-V3i 


2 



cft(A) 



6 6 
5 



-V3i -9+\/3i 
-3+\/3i 2 



s/3i -9-s/3i 
-3-V3i 2 



Operations 

We now briefly illustrate how the eft accelerates and simplifies many operations. Let 
a,§_ € Kfe. Note that 

a + /3 = icft(cft(a) + eft (/})), and 
q o /3 = icf t(cf t(a) cft(/3)). 

In the Fourier space - the output of the eft operation - these operations are both 
0{k) time because they occur between diagonal matrices. Due to the linearity of the 
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eft operation, arbitrary sequences of operations in the Fourier space transform back 
seamlessly, for instance 

(a + P) o (a + P) o . . . o (a + p) = icf t((cf t(a) + cf t(P)) j ). 



j times 

But even more importantly, these simplifications generalize to matrix-based operations 
too. For example, 

iox = icf t(cf t(A) cf t(x)). 
In fact, in the Fourier space, this system is a series of independent matrix vector products: 



cf t(A) cf t(x) 









xi 














Xfc _ 





AfcXfc 



Here, we have again used Aj and Xj to denote the blocks of Fourier coefficients, or 
equivalently, circulant eigenvalues. The rest of the paper frequently uses this convention 
and shorthand where it is clear from context. This formulation takes 

0(mnk log k + nk log k) + 0(kmn) 

eft and icft matvecs 

operations instead of 0(mnk 2 ) using the circ formulation in the previous section. 

More operations are simplified in the Fourier space too. Let cft(a) = diag [ai, & k ]. 
Because the ctj values are the eigenvalues of circ(a), the following functions simplify: 

abs(a) = i eft (diag [ I |d fc |]), 

a = icf t(diag [37, = icf t(cf t(a)*), and 

angle(a) = icft(diag [cvi/|cn|, & k /\a k \]). 

Complex values in the CFT. A small concern with the icft operation is that it may 
produce complex- valued elements in Kfe. It suffices to note that when the output of a 
sequence of circulant operations produces a real- valued circulant, then the output of icft 
is also real-valued. In other words, there is no problem working in Fourier space instead 
of the real-valued circulant space. This fact can be formally verified by first formally 
stating the conditions under whic h icft produ ces real- valued circulants (icf t(D) is real 



Davis! (|l979h . and then checking that the operations 



if and only if F DF = D* , see 
in the Fourier space do not alter this condition. 

3. 2. Properties 

Representations in Fourier space are convenient for illustrating some properties of 
these operations. 

Proposition 2. The matrix circ(angle(a)) is orthogonal. 
Proof. We have 

circ(angle(a))* circ(angle(a)) -H- 

&iai /\&i\ 2 



angle(a) o angle(a) = icft 



a k a k /\a k \ 2 



I. 



□ 
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Additionally, the Fourier space is an easy place to understand spanning sets and bases 
in K™, as the following proposition shows. 

Proposition 3. Let X £ K™ x ™. Then X spans if and only if circ(X) and cft(X) 
have rank km. Also X is a basis if and only if circ(X) and cft(X) are invertible. 

Proof. First note that rank(cf t(X)) = rank(circ(X )) because eft is a similarity trans- 
formation applied to circ. It suffices to show this result for cft(X), then. Now consider 
y X o a: 

cft(y) = cft(X) eft (a); 



yi 






' x 1 












yfc _ 















Thus, if there is a y that is feasible, then all Xj £ c mxn mu st be rank m. Conversely, if 
cft(X) has rank km then each Xj must have rank m, and any y is feasible. The result 
about the basis follows from an analogous argument. □ 



3.3. Inner products, norms, and ordering 

Wc now return to our inner product and norm to elaborate on the positive-definitcness 
and the triangle inequality. In terms of the Fourier transform, 

(x,y) = icft(cft(y)* cft(x)). 

If we write this in terms of the blocks of Fourier coefficients then 

'y**i 



cf t(x)* cf t(y) 



Yk x k 



For y = x, each diagonal term has the form x*Xj > 0. Consequently, we do consider 
this a positive semi-definite inner product because the output circ((x,y)) is a matrix 
with non-negative eigenvalues. This idea motivates the following definition of element 
ordering. 

Definition 4 (Ordering). Let a, f3 £ Kf.. We write 

a < ft when diag(cft(a)) < diag(cf t(f3)) element-wise, and 
a < ft when diag(cft(a)) < diag(cf t(f3)) element-wise. 

We now show that our inner product satisfies the Cauchy-Schwarz inequality: 

abs (x,y) < ||x|| o ||y|| . 

In Fourier space, this fact holds because ly^Xjl < ||xj|| 1 1 y j 1 1 follows from the standard 
Cauchy-Schwarz inequality. Using this inequality, we find that our norm satisfies the 
triangle inequality: 

||x + yjf = (x + y,x + y) < (x,x) +2o ||x|| o ||y|| + (y,y) = (||x|| + ||y||) 2 . 
In this expression, the constant 2 is twice the multiplicative identify, that is 2 = {2 ... 0}. 
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4. Eigenvalues and Eigenvectors 



With the results of the previous fe w sections, we can now state and analyze an 
eigenvalue problem in circulant algebra. Bramanl ( To appear ) investigated these already 
and proposed a decomposition approach to compute them. We offer an extended analysis 
that addresses a few additional aspects. Specifically, we focus on a canonical set of 
eigenpairs. 

Recall that eigenvalues of matrices are the roots of the characteristic polynomial 
det(A - A J) = 0. Now let A e K£ x ™ and A e K k . The eigenvalue problem does not 
change: 

det(A- Ao I) =0. 

(As an aside, note that the standard properties of the determinant hold for any matrix 
over a commutative ring with identity; in particular, the Cayley-Hamilton theorem holds 
in this algebra.) The existence of an eigenvalue implies the existence of a corresponding 
eigenvector x 6 Thus, an eigenvalue and eigenvector pair in this algebra is 

A o x = A o x. 



Just like the matrix case, these eigenvectors can be rescaled by any constant a £ K^: 
Aoaox = Aoaox. In terms of normalization, note that \\P ° x|| = ||x|| if circ(/3) is 
an orthogonal circulant. This follows most easily by noting that 



o x <-> 



\ 1/2 

circ(/3)* circ(xi)* circ(xi) circ(/?) -H> ||x| 



because circulant matrices commute and circ(/3) is orthogonal by construction. For this 
reason, we consider orthogonal circulant matrices the analogues of angles or signs, and 
normalized eigenvectors in the circulant algebra can be rescaled by them. (Recall that 
we showed that angle(a) is an orthogonal circulant in Section [3]) 

The Fourier transform offers a c onvenient decoupling procedure to compute eigenval- 
ues and eigenvectors, as observed bv lBraman ( To appear ). Let A e 
and A be an eigenvalue and eigenvector pair: Aox = xoA and det(A — Ao7) 
it is straightforward to show that the Fourier transforms cft(A), cft(x), and cft(A) 
decouple as follows: 

cf t(A ox) = cf t(x o A); 
cft(A) cft(x) = cft(x) cft(A); 



™ xn and letxeK™ 
0. Then 



A h J 
Aix 





"Ai 








A fc . 



A k i k 



Aiici 



where Xj € X(Aj) and Xj ^ 0. The last equation follows because 

cf t(det(A - Ao/)) = diag [ det(Ai-Aii), dct(A k -\ k i) ] = 0. 
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The decoupling procedure we just described shows that any eigenvalue or eigenvector 
of A must decompose into individual eigenvalues or eigenvectors of the eft-transformed 
problem. This illustrates a fundamental difference from the standard matrix algebra. 
For standard matrices, requiring det(A — XI) = and finding a nonzero solution x for 
Ax = Ax are equivalent. In contrast, the determinant and the eigenvector equations are 
not equivalent in the circulant algebra: 4ox = xoA actually has an infinite number of 
solutions A. For instance, set xi, Ai to be an eigenpair of A\ and Xj = for j > 1, then 
any value for Xj solves i 4ox = xoA. However, only a few of these solutions also satisfy 
det(A- Ao7) = 0. 

Eigenvalues of matrices in K^ xn have some interesting properties. Most notably, a 
matrix may have more than n eigenvalues. As a special case, the diagonal elements of a 
matrix are not necessarily the only eigenvalues. We demonstrate these properties with 
an example. 



Example 5. For the diagonal matrix 



{2 31} {00 0} 
{0 00} {31 1} 



Thus, 



Ax = 



6 
5 



Ao = 



-iVs o 
2 



A, = 



iVs o 
2 



Ai = ic/t(diag [6 22])= (1/3) {10 4 4} 
A 3 = ic/i(diag [ 6 -iV3 tVs ]) = {231} 



A 2 = ic/i(diag [5 -1V3 iVS]) = (1/3) {s 2 2} 
A4 = ic/i(diag [522])= (1/3) {311}. 



The corresponding eigenvectors are 



2?i 



X3 = 



{1/3 1/3 1/3} 
{2/3 -1/3 -1/3} 

{100} 
{0 0} 



X2 



X 4 = 



{2/3 -1/3 -1/3} 
{1/3 1/3 1/3} 

{0 0} 

{10 0} 



There are still more eigenvalues, however. The four eigenvalues above all correspond 
to elements in with real-valued entries. We can combine the eigenvalues of the Aj 's 
to produce complex-valued elements in that are also eigenvalues. These are 



A5 = ic/t(diag [6 -1V3 2]) 
A7 = ic/t(diag [5 -1V3 2]) 



A6 = ic/t(diag [6 2 iVS]) 
A 8 = ic/t(diag [5 2 iVs]). 



For completeness and further clarity, let us extend this example a bit by presenting also 



the eigenvalues of the non-diagonal matrix from Example^ Let A = 
The eft produces: 



{2 3 1} {8 -2 0} 
{-2 02} {31 1} 



6 6 
5 



V3 
-3 + iV3 



-9 + W3 
2 



A, = 



iV3 -9 - iV3 
-3 + iV3 2 
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The numerical eigenvalues of A\ are {6,5}; of Ai are {-0.0899 + 6. 4282z, 2.0899 — 
4.6962i}; and of A 3 are {-0.0899 - 6.4282«, 2.0899 + 4.6962?}. The real-valued eigen- 
values of A. are 

Ai = {l.9401 -1.6814 5.7413} A2 = {3.0599 3.6814 -1.7413} 
A3 = {3.3933 4.0147 -1.4080} A4 = {1.6067 -2.0147 5.4080} . 

The complex-valued eigenvalues of A. are 

A 5 = {4.6966-1.5654i, -0. 7040+1. 9114i, 2. 0073-0. 3461i} 
A 6 = {3.6367+2.1427i, 3. 0373+0. 3980i, -0. 6740-2. 5407i} 
A7 = {4.3633-1.5654i, -1. 0373+1. 9114i, 1. 6740-0. 3461i} 
A 8 = {3.3034+2.14274, 2. 7040+0. 3980i, -1. 0073-2. 5407i} . 

We now count the number of unique eigenvalues and eigenvectors, using the decou- 
pling procedure in the Fourier space. To simplify the discussion, let us only consider the 
case where each Aj has simple eigenvalues. Consider an A. £ K? x " with this property, 
and let m,j be the number of unique eigenvalues and eigenvectors of Aj . Then the number 
of unique eigenvalues of A is given by the number of unique solutions to det( A-XoI) = 
which is rij=i m j- The number of unique eigenvectors (up to normalization) is given by 
the number of unique solutions toAox = xoA, which is also 0.7=1 m j- 

This result shows there are at most n k eigenvalues if A £ Kk is allowed to be complex- 
valued, even when A. £ Kk is real-valued. If A. £ Kk is real-valued, then there are at 
most n^ k+1 " 2 > "real" eigenvalues. For this result, note that icf t(diag [ °i ■■■ a k ]) is real- 
valued if and only if diag [m - at] — F diag [ Q i ■■■ a k }F 2 (|Davisl . ll979h . where F is the 



Fourier transform matrix. This implies a\ is real-valued, and aj = ak-j+i- Applying 
this restriction reduces the feasible combinations of eigenvalues to n^ k+1 ^ 2 ^ . 

Given that there are so many eigenvalues and vectors, are all of them necessary to 
describe A.? We now show this is not the case by making a few definitions to clarify the 
discussion. 

Definition 6. Let A. £ K^ x ™. A canonical set of eigenvalues and eigenvectors is a set 
of minimum size, ordered such that abs(Ai) > abs(A2) > . . . > abs(Afc), which contains 
the information to reproduce any eigenvalue or eigenvector of A. 

In the diagonal matrix from Example[5l the sets {(Ai,xi), (A2,X2)}, {(A3,X3), (A4,X4)}, 
and {(Ai,xi), (A3,X3), (X^^i)} contain all the information to reproduce any eigenpair, 
whereas the set {(Ai,xi), (A3,X3)} does not (it does not contain the eigenvalue 5 of A\). 
In this case, the only canonical set is {(Ai, Xi), (A2, X2)}. This occurs because, by a 
simple counting argument, a canonical set must have at least two eigenvalues, thus the 
set is of minimum size. The choice of Ai and A2 is given by the ordering condition. 
Among all the size 2 sets with all the information, this is the only one with the property 
that abs(Ai) > abs(A2). 

Theorem 7 (Unique Canonical Decomposition). Let A. £ K^ xn where each Aj in the 
cft(A) matrix has distinct eigenvalues with distinct magnitudes. Then A. has a unique 
canonical set of n eigenvalues and eigenvectors. This canonical set corresponds to a basis 
of n eigenvectors, yielding an eigendecomposition 

A = XoAor'. 
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Proof. Because all of the eigenvalues of each Aj are distinct, with distinct magnitudes, 
there are nk distinct numbers. This implies that any canonical set must have at least n 
eigenvalues. 

Let A^ } be the ith eigenvalue of Aj ordered such that \\f ) \ > \xf ] \ > ... > \X { " } \. 
Then Aj = icft(diag X± , . . . , X^ ) is a canonical set of eigenvalues. We now show 

that this set constitutes an eigenbasis. Let Aj = XjAjJtj be the eigendecomposition 
using the magnitude ordering above. Then X = icf t(diag [aj, . . . , Xk] ) and A = 
icf t(diag [Afc, . . . , Afe] ) is an eigenbasis because the matrix X satisfies the properties 
of a basis from Theorem [3J Note that Aj ( j = Aj. 

Finally, we show that the set is unique. In any canonical set abs(Ai) > abs(Aj) for 
i > 1. In the Fourier space, this implies |A^| > |A^ |. Because all of the values 

are unique, there is no choice for A,) in a canonical set and we have |a5 1} | > |Aj°|,» > 1. 
Consequently, Ai is unique. Repeating this argument on the remaining choices for Aj 
shows that the entire set is unique. □ 

Remark 1. If Aj has distinct eigenvalues but they do not have distinct magnitudes, 
then A. has an eigenbasis but the canonical set may not be unique, because Aj may have 
two distinct eigenvalues with the same magnitude. 

Next, we show that the eigendecomposition is real-valued under a surprisingly mild 
condition. 

Theorem 8. Let A. £ K^ xn be real-valued with diagonalizable Aj matrices. If k is odd, 
then the eigendecomposition X o Ao X 1 is real-valued if and only if Ai has real-valued 
eigenvalues. If k is even, then X o A o X 1 is real-valued if and only if Ai and ^4.^/2+1 
have real-valued eigenvalues. 

Proof. First, if A. has a real- valued eigendecomposition, then we have that X\ is real 
and also that X^/2+i is real when k is even. Likewise, Ai is real and Afc/ 2 +i is real when 
k is even. Thus, A\ (and also A^/2+i when k is even) have real-valued eigenvalues and 
vectors. 

When Ai (and A k / 2 +i for k even) have real-valued eigenvalues and vectors, then note 
that we can choose eigenvalues and eigenvectors of the other matrices Aj, which may 
be complex, in complex-conjugate pairs so as to satisfy the condition for a real-valued 
inverse Fourier transforms. This happens because when A. is real, then A\ is real and 



Aj = Ak-j+2 by the properties of the Fourier transform ([Davis, Il979h . Thus for each 

eigenpair Aj,Xj of Aj, the pair Aj,Xj is an eigenpair for Ak-j+2- Consequently, if we 
always choose these complex conjugate pairs for all j besides j = 1 (and j = k/2 + 1 for 
k even), then the result of the inverse Fourier transform will be real- valued. □ 

Finally, we note that if the scalars of a matrix are padded with zeros to transform 
them into the circulant algebra, then the canonical set of eigenvalues are nothing but 
tuples that consist of the eigenvalues of the original matrix in the first entry, padded with 
zeros as well. To justify this observation, let A. g K^ xn have Aij — {Gij, 0, 0} for a 
matrix G € M nx ". Also, let Ai, . . . , A m (m < n) be the eigenvalues of G ordered such 
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Require: A, x^,r 
1: {Kept for alignment} 



2: X 



(0) 



(0) 



A0)\ 



3: for k = 1 , . . . until convergence do 
5: <- ||y( fe ) || 



if converged then 

return x( fc ) 
end if 
end for 



Require: A,xW,t 
<- eft (J 

(o) . -.(0) 



l: A <- cft(A),X (0) 



2: X 



x (0) x (0) 



cft(x(°)) 

1/2 



3: 
4: 

5: 

(5: 
7: 



for A; = 1 , . . . until convergence do 



- (*0 



R 



(k) 



AX {k ~ 1] 

~ <k) -.(fe)- 1 / 2 



if converged then 



- (fe) 

return icft(_X" ) 
end if 
end for 



Figure 3: The power method in the circulant algebra (left) and the power method in the circulant 
algebra after transformation with the fast Fourier transform (right). We address convergence criteria in 
Section 15.11 



that |Ai| > |A 2 | > • • ■ > |A m |. Then cft^A^) = diag[G»,i, G i:j } and thus Aj = G for 
all j. Thus, we only need to combine the same m eigenvalues of each Aj to construct 
eigenvalues of A. For the eigenvalues Ai, we have cft(A) = diag[Ai, a 4 ], thus the 
given set is canonical because of the same argument used in the proof of Theorem [7J 

We end this section by noting that much of the above analysis can be generalized to 
non-simple eigenvalues and vectors using the Jordan canonical form of the Aj matrices. 



5. The power method and the Arnoldi method 

In what follows, we show that the power method in the circulant algebra computes 
the eigenvalue Ai in the canonical set of eigenvalues. This result shows how the circulant 
algebra matches the behavior of the standard power method. As part of our analysis, we 
show that the power method decouples into k independent power iterations in Fourier 
space and is equivalent to a subspace iteration method. Second, we demonstrate the 
Arnoldi method in the circulant algebra. In Fourier space, the Arnoldi method is also 
equivalent to the Arnoldi algorithm on independent problems, and it also corresponds to 
a particular block Arnoldi procedure. 

5.1. The power method 

Please see the left half of Figure [3] for the sequence of operations in the power method 
in the circulant algebra. In fact, it is not too different from the standard power method 
in Figure Q] We replace Ax with Aox and use the norm and inverse from Section [3] 
We'll return to the convergence criteria shortly. As we show next, the algorithm runs k 
independent power methods in Fourier space. Thus, the right half of Figure |3] shows the 
equivalent operations in Fourier space. 
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To analyze the power method, consider the key iterative operation in the power 
method when transformed into Fourier space: 



cft(Aoxo (||A.ox||) 1 ) 

= eft (A) cft(x)(cft(x)* cft(x))" 1/2 



Aixi 



Aixi 



A fc i fc 



-1/2 



Now, 



Aiii 



A k i 



Aiii 



-1/2 



ij A 1 Aiii 



Aiii 



£A fc A fc i fc 



-1/2 



Thus 



cftQAoxo (||Aox||) 1 ) = 



Aiii/ Aiii 



A fc i fc / Aiii 



The key iterative operation, Aoxo (||A o x||) , corresponds to one s£ep of the standard 
power method on each matrix Aj . From this derivation, we arrive at the following theo- 
rem, whose proof follows immediately from the convergence proof of the power method 
for a matrix. 

Theorem 9. Let A £ K^ xn have a canonical set of eigenvalues X±, . . . , X n where |AjJ > 
|A2 1 j then the power method in the circulant algebra convergences to an eigenvector xj 
with eigenvalue Ai . 

A bit tangcntially, a n eigenpair in the Fourier sp ace is a simple instance of a multivari- 
ate eigenvalue problem ( Chu and Wattersonl . Il993) . The general multivariate eigenvalue 



problem i s A^x, = AjXj i = 1, . . ., whereas we study the same system, albeit 
diagonal. IChu and Wattersonl |l993) did study a power method for the more general 
problem and showed local convergence; however our diagonal situation is sufficiently 
simple for us to state stronger results. 

Convergence Criteria. A simple measure such as \\x^ — x^" -1 ) || < r, with t = {tO ... o} 
will not detect convergence. As mentioned in the description of the standard power 
method in Figure [TJ this test can fail when the eigenvector changes angle. Here, we have 
the more general notion of an angle for each element, and eigenvectors are unique up 
to a choice of angle. Thus, we first normalize angles before comparing the we use the 
convergence criteria 



angle(xi 



(k) 



^oxW 



angle(x^ ) ox 



,0-1) 



< T. 



(0) 
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In the Fourier space, this choice requires that all of the independent problems have con- 
verged to a tolerance of r, which is a viable practical choice. An alternative convergence 
criteria is to terminate when the eigenvalue stops changing, although this may occur 
significantly before the eigenvector has converged. 

Subspace iteration. We now show that the power method is equivalent to subspace iter- 
ation in Fourier space. Subspace iteration is also known as "orthogonal iteration" or the 
"block-power method." Given a starting block of vectors X^°\ the iteration is 

Y^AX {k) , x {k+1 \R {k+1] =qr(T). 

On the surface, there is nothing to relate this iteration to our power method, even 
in Fourier space. The relationship, however, follows because all of our operations in 
Fourier space occur with block- diagonal matrices. Note that for a block-diagonal matrix 

- (k) 

of vectors, which is what X is, the QR factorization just normalizes each column. In 
other words, the result is a diagonal matrix R. This simplification shows that steps[5]|n]in 
the Fourier space algorithm are equivalent to the QR factorization in subspace iteration. 

Breakdown. One problem with this iterative approach is that it can encounter "zero 
divisors" as scalars when running these algorithms. These occur when the matrices in 
Fourier space are not invertible. We have not explicitly addressed this situation and note 
that the same issues arise in block methods when some of the quantities become singular. 
The analogy with the block method may provide an appropriate solution. For example, 

if the scalar gS k ^ is a zero-divisor, then we could use the QR factorization of '" - as 
suggested by the equivalence with subspace iteration - instead. 

5.2. The Arnoldi process 

The Arnoldi method is a cornerstone of modern matrix computations. Let A be an 
n x n matrix with real valued entries. Then the Arnoldi method is a technique to build 
an orthogonal basis for the Krylov subspace 

JC t {A, v) = spanjv, Av, . . . , A t_1 v}, 

where v is an initial vector. Instead of using this power basis, the Arnoldi process 
derives a set of orthogonal vectors that span the same space when computed with exact 
arithmetic. The standard method is presented in Figure SJa). From this procedure, we 
have the Arnoldi decomposition of a matrix: 

AQ t — Q t+1 H t+1}t 

where Q t is an n x t matrix, and Ht+i.t is a (i+1) x t upper Hessenberg matrix. Arnoldi's 
orthogon al subspaces Q enable efficient algorithms for both solvin g large scale li near 
systems ( Krvlovl Il931l ) and computing eigenvalues and eigenvectors ( Arnoldi Il951 ) . 

Using our repertoire of operations, the Arnoldi method in the circulant algebra is 
presented in Figure BJb) . The circulant Arnoldi process decoupled via the eft is also 
shown in Figure Etc). 

We make three observations here. First, the decoupled (eft) circulant Arnoldi process 
is equivalent to individual Arnoldi processes on each matrix Aj . This follows by a similar 
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(a) Arnoldi for R nx " 
Require: A, b, t 

I: 

2: 

3: qi <-b/||b|| 

4: for j = 1, . . . ,t do 

5: Z <- Aqi 

6: for i = 1, . . . , i do 

7: //,., <- q*z 

8: 

9: Z <- Z - ffij-q, 

10: end for 

11: JT i+1J <- ||z|| 
12: 

13: q J + i <- z/Hj+ij 
14: 

15: end for 

Figure 4: Arnoldi methods. Algorithm (a) shows the standard Arnoldi process. Algorithm (b) shows 
the Arnoldi process in the circulant algebra, and Algorithm (c) shows the set of operations in (b) but 
expressed in the Fourier space. 

analysis used to show the decoupling result about the power method. The verification of 
this fact for the Arnoldi iteration is a bit more tedious and thus we omit this analysis. 
Second, the same decoupled process is equivalent to a block Arnoldi process. This also 
follows for the same reason the equivalent result held for the power method: the QR 
factorization of a block-diagonal matrix-of- vectors is just a normalization of each vector. 
Third, we produce an Arnoldi factorization: 

AoQ t = Q t+1 off t +i,(. 

In fact, this outcome is a corollary of the first property and follows from applying icf t 
to the same analysis. 

This discussion raises an interesting question, why iterate on all problems simultane- 
ously? One case where this is advantageous is with sparse problems; and we return to 
this issue in the concluding discussion (Section [8]). 

6. A Matlab package 

The Matlab environment is a convenient playground for algorithms involving ma- 
trices. We have extended it with a new class implementing the circulant algebra as a 
native Matlab object. The name of the resulting package and class is camat: circulant 
algebra matrix. While we will show some non-trivial examples of our package later, let 
us start with a small example to give the flavor of how it works. 

A = cazeros (2 , 2 , 3) ; °L creates a camat type 



(b) Arnoldi for K" xn 
Require: A,h,t 

1: 
2: 

3: qi <- b o 1 1 b 1 1 1 

4: for j = 1, . . . ,t do 

5: z <— A. o q j 

6: for i = 1, . . . , j do 

7: H hJ 4- (q l; z) 
8: 

9: Z*-Z- Hij o q, 

10: end for 

11: Hj +1J <- ||z|| 
12: 

13: Cjj+i <- Z o //,.',, 
14: 

15: end for 



(c) Unrolled Arnoldi for K£ xn 
Require: A,h,t 

1: A <- cf t(A) 

2: B <- Cf t(b) 

3: Q 1 <- B(B*B)-V 2 

4: for j = 1, . . . ,t do 

5: Z <- AQj 

6: for i = 1, . . . , j do 

7: Hij 4- Q*Z 

8: H, , 4- icf t; 77,.,) 

9: Z 4r- Z Q.Hi j 

10: end for 

11: H j+1J <- {ZZfll 

12: !/, I , <- icft(H j+1<j ) 

13: Q. 4- Z1I t 

14: q j+1 <- icft(Q J + 1 ) 

15: end for 
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A(l,l) = cascalar( [2,3, 1] ) ; A(l,2) = cascalar( [8,-2,0] ) ; 
A(2,l) = cascalar( [-2,0,2] ) ; A(2,2) = cascalar ( [3, 1 , 1] ) ; 
eig(A) % compute eigenvalues as in Example 2; 

The output, which matches the non-diagonal matrix in Example [3 is: 

ans = 

(:,:,1) = % the first eigenvalue 
1.9401 
-1.6814 
5.7413 

(:,:,2) = 7, the second eigenvalue 
3.0599 
3.6814 
-1.7413 

Internally, each element A £ K™ xn is stored as a k x n x m array along with its eft 
transformed data. Each scalar is stored by the k parameters defining it. To describe this 
storage, let us introduce the notation 

= circ(a)ei, 

to label the vector of k parameters explicitly. Thus, we store vec(a) for a G Kj,. This 
storage corresponds to storing each scalar consecutively in memory. The matrix is 
then stored by rows. We store the data for the diagonal elements of the eft transformed 
version in the same manner; that is, diag(cft(a)) is stored as k consecutive complex- 
valued scalars. The organization of matrices and vectors for the eft data is also by row. 
The reason we store the data by row is so we can take advantage of Matlab's standard 
display operations. 

At the moment, our implementation stores the elements in both the standard and 
Fourier transformed space. The rationale behind this choice was to make it easy to 
investigate the results in this manuscript. Due to the simplicity of the operations in 
the Fourier space, most of the functions on camat objections use the Fourier coefficients 
to compute a result efficiently and then compute the inverse Fourier transform for the 
vec representation. Hence, rather than incurring for the Fourier transform and inverse 
Fourier transform cost for each operation, we only incur the cost of the inverse transform. 
Because so few operations are easier in the standard space, we hope to eliminate the 
standard vec storage in a future version of the code to accelerate it even further. 

We now show how the overloaded operation eig works in Figure [5] This procedure, 
inspired by Theorem [8j implements the process to get real- valued canonical eigenvalues 

and eigenvectors of a real- valued matrix in the circulant algebra. The slice Af ( j , : , : ) is 

-t 

the matrix A ■ . Here, the real-valued transpose results from the storage-by-rows instead 
of the storage-by-columns. The code proceeds by computing the eigendecomposition of 
each Aj with a special sort applied to produce the canonical eigenvalues. After all of the 
eigendecompositions are finished, we need to transpose their output. Then it feeds them 
to the if ft function to generate the data in vec form. 

In a similar manner, we overloaded the standard assignment and indexing opera- 
tions e.g. a = A(i,j); A (1,1) = a; the standard Matlab arithmetic operations +, -, 
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function [V,D] = eig(A) 

7» CAEIG The eigenvalue routine in the circulant algebra 

Af = A.fft; k = size(Af,l); 7« extract data from object 

if any(imag(A.data(: ))) , error (' specialized for real values'); end 

[Vf,Df] = deal (zeros (Af )) ; °L allocate data of size (k,n,n) 

[Vf (1, : , :) ,Df (1, : , :)] = sortedeig(squeeze(Af (1 ,:,:)).') ; 

for j=2:floor(k/2)+l 

[Vf (j , : , :) ,Df (j , : , :)] = sortedeig(sqeeze(Af (j ,:,:)). ') ; 

if j~=k/2+l 7 skip last when k is even 

Vf (k-j+2, : , :) = conj (Vf (j ,:,:)) ; Df (k-j+2 , : , : ) = conj (Df (j ,:,:)) ; 

end 
end 

% transpose all the data back. 

for j=l:k, Vf(j,:,:) = Vf(j, :,:).'; Df(j,:,:) = Df(j, :,:).'; end 

V = camatcft(if ft(Vf ) ,Vf ) ; % create classed output 
D = camatcft(if ft(Df ) ,Df ) ; 

function [V,D] =sortedeig(A) 

[V,D] = eig(A); d = diag(D) ; [ignore p] = sort (-abs (d) ) ; 

V = V(:,p); D = D(p,p); 7. apply the sort 

Figure 5: The implementation of the eigenvalue computation in our package. Please see the discussion 
in the text. 

for iter=l :maxiter 
Ax = A*x; 
lambda = x'*Ax; 
x2 = (1./ norm(Ax))*Ax; 

delta = mag(norm(l./angle(x(l))*x-l./angle(x2(l))*x2)) ; 
if delta<tol, break, end 
end 

Figure 6: The implementation of the power method using our package. 

*, /, \; and the functions abs , angle, norm, conj, diag, eig, hess, mag, norm, 
numel, qr, rank, size, sqrt , svd. 

All of these operations have been mentioned or are self explanatory, except mag. It 
is a magnitude function, and we discuss it in detail in |Appcndix A| 

Using these overloaded operations, implementing the power method is straightfor- 
ward; see Figure [6] We note that the power method and Arnoldi methods can be further 
optimized by implementing them directly in Fourier space. This remains as an item for 
future work. 

7. Numerical examples 

In this section, we present a numerical example using the code we described in Sec- 
tion [5J The problem we consider is the Poisson equation on a regular grid with a mixture 
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of periodic and fixed boundary conditions: 

-Au(x,y) = f(x,y) u(x,0) = u(x, l),u(0,y) = y(l,y) = (x,y) e [0,1] x [0,1]. 

Consider a uniform mesh and the standard 5-point discrete Laplacian: 

-Au(xi,yj) fts -u{xi-uVj) - uixi^j-!) + Au(x l ,y j ) - u(x i+1 ,yj) - u(xi,y j+1 ). 

After applying the boundary conditions and organizing the unknowns of u in y-major 
order, an approximate solution u is given by solving an N(N — 1) x N(N — 1) block- 
tridiagonal, circulant-block system: 

C -I 
-I c 

. . J 

-I c 





u(xi,-) 








" 4 -1 -1" 




U(X 2 , ■) 






c = 


-1 4 '•• 




Vi{x N -\, ■)_ 




_f(x N -i, •)_ 




'•• '•■ -1 








-1 -1 4_ 


u 




f 




NxN 



that is, Au = f . Because of the circulant-block structure, this system is equivalent to 

iou = f 

where A. is a,n N — lx N — 1 matrix of Kjv elements, u and f have compatible sizes, and 

A = circ(A) u = vec(u) f = vec(f). 
We now investigate this matrix and linear system with N = 50. 

7.1. The power method 

We first study the behavior of the power method on A. The canonical eigenvalues of 
A are 

Xj = {4+2 cosO'tt/AT), -1,0,. ..,0,-1} . 

To sec this result, let A(/i) = {^,-1,0,.. .,0,-1} . Then 

"(4-M)°l 



(A-X(fx)ol) 



-lol (4-/i)ol 



-lol 
lol (4-/i)ol 



The canonical eigenvalues of A. — A(/i) o / can be determined by choosing /1 to be an 
eigenvalue of T = tridiag(— 1, 4, —1). These are given by setting /1 = 4 + 2cos(j7r /N), 
where each choice j — 1, . . . , N— 1 produces a canonical eigenvalue Xj. From these canon- 
ical eigenvalues, we can estimate the convergence behavior of the power method. Recall 
that the algorithm runs independent power methods in Fourier space. Consequently, 
these rates are given by A2/A1 for each matrix Aj. To state these ratios compactly, let 
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71 = 4 + 2cos(tt/7V) and 72 = 4 + 2 cos(27r/iV); also let Sj = 2cos(-?r + 27r(j - 1)/N). 
For N even, 

cft(Ai) = diag [ti+<5i, — >ti+*jv] 
cft(Aa) = diag [72+61; •••,72+^ ] 

Thus, the convergence ratio for Aj is (72 + ^)/(7i + The largest ratio (fastest 
converging) corresponds to the smallest value of Sj, which is Si. The smallest ratio 
(slowest converging) corresponds to the largest value of Sj, which is S N / 2 +i in this case. 
(This choice will slightly change in an obvious manner if N is odd.) Evaluating these 
ratios yields 

min b&l = Ii±ii = 2 + 2cos(2 7 r/iV) 
i Ai(lj) 7i+^i 2 + 2cos(7r/iV) 

\ 2 {Aj) 72+^/2+1 6 + 2cos(27r/iV) 

max 7T— = = — — 7 — tttt slowest). 

i Xi(Aj) li+S N/2+1 6 + 2cos(7r/iV) 

Based on this analysis, we expect the eigenvector to converge linearly with the rate 
6 6+2°os(tv/n) • the standard theory for the power method, expect the eigenvalues to 
converge twice as fast. 

Let p be the eigenvector change measure from equation ([5]). In Figure [71 we first show 
how the maximum absolute value of the Fourier coefficients in p behaves (the red line). 
Formally, this measure is ||cft(p)|L, i.e., the maximum element in the diagonal matrix. 
We also show how each Fourier component of the eigenvalue converges to the Fourier 
components of Ai (each gray line). That is, let //W be the Rayleigh quotient xW oA.oxS' 1 ' 
at the ith iteration. Then these lines are the N values of diag(cf t(abs(/xW _^ : ))), The 
results validate the theoretical predictions, and the eigenvalue does indeed converge to 
Ai. 



7.2. The Arnoldi method 

We next investigate computing u using the Arnoldi method applied to A.. In this 
case, f(x,y) to be 1 at £25,2/2 and elsewhere. This corresponds to a single non-zero in 
vec(f) with value 1/N 2 . With this right-hand side, the procedure we use is identical to 
an unoptimized GMRES procedure. Given a t-step Arnoldi factorization starting from 
f , we estimate 

u (t) ~Qt° argmin ||i?t+i,t y - P §i|| , 

yeK fc 

where (3 = ||f || . We solve the least-squares problem by solving each problem inde- 
pendently in the Fourier space - as has become standard throughout this paper. Let 
p = llf — A o uw|| . Figure [S] shows (in red) the magnitude of the residual as a func- 
tion of the Arnoldi factorization length t, which is ||cft(/j)|| . The figure also shows (in 
gray) the magnitude of the error in the jth Fourier coefficient; these lines are the N 
values of diag(cf t(||u — u^' || )). In Fourier space, these values measure the error in each 
individual Arnoldi process. 

What the figure shows is that th e residual su ddenly converges at the 26th iteration. 
This is in fact theoretically expected (|Saadl . 120031 ). because each matrix Aj has A^/2 + 1 = 
26 distinct eigenvalues. In terms of measure the individual errors (the gray lines), some 
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Figure 7: The convergence behavior of the power 
method in the circulant algebra. The gray lines 
show the error in the each eigenvalue in Fourier 
space. These curves track the predictions made 
based on the eigenvalues as discussed in the text. 
The red line shows the magnitude of the change 
in the eigenvector. We use this as the stopping 
criteria. It also decays as predicted by the ra- 
tio of eigenvalues. The blue fit lines have been 
visually adjusted to match the behavior in the 
convergence tail. 



Figure 8: The convergence behavior of a GM- 
RES procedure using the circulant Arnoldi pro- 
cess. The gray lines show the error in each 
Fourier component and the red line shows the 
magnitude of the residual. We observe poor 
convergence in one Fourier component; until the 
Arnoldi basis captures all of the eigenvalues af- 
ter N/2 + 1 = 26 iterations. These results show 
how the two computations are performing indi- 
vidual power methods or Arnoldi processes in 
Fourier space. 



converge rapidly, and some do not seem to converge at all until the Arnoldi process 
completes at iteration 26. This exemplifies how the overall behavior is governed by the 
worst behavior in any of the independent Arnoldi processes. 



8. Summary 



We have extended the circulant algebra, introduced bv lKilmer et al. (2008). with new 



operations to pave the way for iterative algorithms, such as the power method and the 
Arnoldi iteration that we introduced. These operations provided key tools to build a 
Matlab package to investigate these iterative algorithms for this paper. Furthermore, 
we used the fast Fourier transform to accelerate these operations, and as a key analysis 
tool for eigenvalues and eigenvectors. In the Fourier space the operations and algorithms 
decouple into individual problems. We observed this for the inner product, eigenvalues, 
eigenvectors, the power method, and the Arnoldi iteration. We also found that this 
decoupling explained the behavior in a numerical example. 

Given that decoupling is such a powerful computational and analytical tool, a natural 
question that arises is when it is useful to employ the original circulant formalism, rather 
than work in the Fourier space. For dense computations, it is likely that working entirely 
in Fourier space is a superior approach. However, for sparse computations, such as the 
system Aou = f explored in Section [3 such a conclusion is unwarranted. That example 
is sparse both in the matrix over circulants, and in the individual circulant arrays. When 
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thought of as a cube of data, it is sparse in any way of slicing it into a matrix. After this 
matrix A is transformed to the Fourier space, it loses its sparsity in the third-dimension; 
each sparse scalar Aij becomes a dense array. In this case, retaining the coupled nature of 
the operations and even avoiding most of the Fourier domain may allow better scalability 
in terms of total memory usage. 

An interesting topic for future work is exploring other rings besides the ring of circu- 
lants. One obvious candidate is the ring of symmetric circulant matrices. In this ring, 
the Fourier coefficients are always real-valued. Using this ring avoids the algebraic and 
computational complexity associated with complex values in the Fourier transforms. 

We have made all of code and experiments available to use and reproduce our results: 
http : //Stanford. edu/~dgleich/publicat ions/20 11/codes/camat 
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Appendix A. The circulant scalar magnitude 

This section describes another operation we extended to the circulant algebra. Even- 
tually, we replaced it with our ordering (Definition [4]), which is more powerful as we 
justify below. However, it plays a role in our Matlab package, and thus we describe the 
rationale for our choice of magnitude function here. 

For scalars in R, the magnitude is often called the absolute value. Let a, (3 E M. 
The absolute value has the the property \aj3\ = \a\ \/3\. We have already introduced an 
absolute value function, however. Here, we wish to define a notion of magnitude that 
produces a scalar in R to indicate the size of an element. Such a function will have 
norm-like flavor because it must represent the aggregate magnitude of k values with a 
single real-valued number. Thus, finding a function to satisfy la o /? = |a| |/?| exactly is 
not possible. Instead, we seek a function g : h-> R such that 

1- 9(q) = if and only if a = Q, 

2- g(a°£)<g(a)g(0), 

3- g(a + /3) <g(a)+g(0. 

The following result shows that there is a large class of such magnitude functions. 
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Result 1. Any sub- multiplicative matrix norm || A\\ defines a magnitude function g{q) = 
||circ(a)|| . 

This result follows because the properties of the function g are identical to the re- 
quirements of a sub-multiplicative matrix norm applied to circ(a). Any matrix norm 
induced by a vector norm is sub-multiplicative. In particular, the matrix 1, 2, and oo 
norms are all sub-multiplicative. Note that for circulant matrices both the matrix 1 
and oo norms are equal to the 1-norm of any row or column, i.e., ||vec(a)|| 1 is a valid 
magnitude. Surprisingly, the 2-norm of the vector of parameters, that is ||vec(a)|| 2 , 
is not. For a counterexample, let a = {i 2} , /3 = {24}. Then a o (5 = {s 10} and 
||vec(ao/3)|| = V164 > ||vec(a)|| 2 ||vec(/3)|| = >/l00. For many practical computa- 
tions, we use the matrix 2-norm of circ(a) as the magnitude function. Thus, 

\a\ = ||circ(a)|| 2 = Hcftfe)^ . 

This choice has the following relationship with our ordering: 

abs(a) < abs(^) \g\ < . 

We implement this operation as the mag function in our Matlab package. 
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